home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Language/OS - Multiplatform Resource Library
/
LANGUAGE OS.iso
/
cpp_libs
/
lin_alg.lha
/
lin_alg
/
vmatrix.cc
< prev
next >
Wrap
C/C++ Source or Header
|
1993-08-08
|
14KB
|
460 lines
// This may look like C code, but it is really -*- C++ -*-
/*
************************************************************************
*
* Verify Operations on Matrices
*
************************************************************************
*/
#include "LinAlg.h"
#include <math.h>
#include <builtin.h>
#include <ostream.h>
/*
*------------------------------------------------------------------------
* Service validation functions
*/
static void verify_element_value(const Matrix& m,const REAL val)
{
register imax = 0, jmax = 0;
register double max_dev = 0;
register int i,j;
for(i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
{
register double dev = abs(m(i,j)-val);
if( dev >= max_dev )
imax = i, jmax = j, max_dev = dev;
}
if( max_dev == 0 )
return;
else if( max_dev < 1e-5 )
message("Element (%d,%d) with value %g differs the most from what\n"
"was expected, %g, though the deviation %g is small\n",
imax,jmax,m(imax,jmax),val,max_dev);
else
_error("A significant difference from the expected value %g\n"
"encountered for element (%d,%d) with value %g",
val,imax,jmax,m(imax,jmax));
}
static void verify_matrix_identity(const Matrix& m1, const Matrix& m2)
{
register imax = 0, jmax = 0;
register double max_dev = 0;
register int i,j;
are_compatible(m1,m2);
for(i=m1.q_row_lwb(); i<=m1.q_row_upb(); i++)
for(j=m1.q_col_lwb(); j<=m1.q_col_upb(); j++)
{
register double dev = abs(m1(i,j)-m2(i,j));
if( dev >= max_dev )
imax = i, jmax = j, max_dev = dev;
}
if( max_dev == 0 )
return;
if( max_dev < 1e-5 )
message("Two (%d,%d) elements of matrices with values %g and %g\n"
"differ the most, though the deviation %g is small\n",
imax,jmax,m1(imax,jmax),m2(imax,jmax),max_dev);
else
_error("A significant difference between the matrices encountered\n"
"at (%d,%d) element, with values %g and %g",
imax,jmax,m1(imax,jmax),m2(imax,jmax));
}
/*
*------------------------------------------------------------------------
* Test allocation functions and compatibility check
*/
static void test_allocation(void)
{
cout << "\n\n---> Test allocation and compatibility check\n";
Matrix m1(4,20); m1.set_name("Matrix m1");
Matrix m2(1,4,1,20); m2.set_name("Matrix m2");
Matrix m3(1,4,0,19); m3.set_name("Matrix m3");
Matrix m4(m1); m4.set_name("Matrix m4");
cout << "The following matrices have been allocated\n";
m1.info(), m2.info(), m3.info(), m4.info();
cout << "\nStatus information reported for matrix m3:\n";
cout << " Row lower bound ... " << m3.q_row_lwb() << "\n";
cout << " Row upper bound ... " << m3.q_row_upb() << "\n";
cout << " Col lower bound ... " << m3.q_col_lwb() << "\n";
cout << " Col upper bound ... " << m3.q_col_upb() << "\n";
cout << " No. rows ..........." << m3.q_nrows() << "\n";
cout << " No. cols ..........." << m3.q_ncols() << "\n";
cout << " No. of elements ...." << m3.q_no_elems() << "\n";
cout << " Name " << m3.q_name() << "\n";
cout << "\nCheck matrices 1 & 2 for compatibility\n";
are_compatible(m1,m2);
cout << "Check matrices 1 & 4 for compatibility\n";
are_compatible(m1,m4);
#if 0
cout << "Check matrices 1 & 3 for compatibility\n";
are_compatible(m1,m3);
#endif
cout << "m2 has to be compatible with m3 after resizing to m3\n";
m2.resize_to(m3);
are_compatible(m2,m3);
Matrix m5(m1.q_nrows()+1,m1.q_ncols()+5); m5.set_name("Matrix m5");
cout << "m1 has to be compatible with m5 after resizing to m5\n";
m5.info();
m1.resize_to(m5.q_nrows(),m5.q_ncols());
are_compatible(m1,m5);
}
/*
*------------------------------------------------------------------------
* Test uniform element operations
*/
static void test_element_op(const int rsize, const int csize)
{
const double pattern = 8.625;
register int i,j;
cout << "\n\n---> Test operations that treat each element uniformly\n";
Matrix m(-1,rsize-2,1,csize);
Matrix m1(m);
cout << "Writing zeros to m...\n";
for(i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
m(i,j) = 0;
verify_element_value(m,0);
cout << "Clearing m1 ...\n";
m1.clear();
verify_element_value(m1,0);
cout << "Comparing m1 with 0 ...\n";
assure(m1 == 0, "m1 is not zero!");
cout << "Writing a pattern " << pattern << " by assigning to m(i,j)...\n";
for(i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
m(i,j) = pattern;
verify_element_value(m,pattern);
cout << "Writing the pattern by assigning to m1 as a whole ...\n";
m1 = pattern;
verify_element_value(m1,pattern);
cout << "Comparing m and m1 ...\n";
assure(m == m1, "m and m1 are unexpectedly different!");
cout << "Comparing (m=0) and m1 ...\n";
assert(!(m.clear() == m1));
cout << "\nClear m and add the pattern\n";
m.clear();
m += pattern;
verify_element_value(m,pattern);
cout << " add the doubled pattern with the negative sign\n";
m += -2*pattern;
verify_element_value(m,-pattern);
cout << " subtract the trippled pattern with the negative sign\n";
m -= -3*pattern;
verify_element_value(m,2*pattern);
cout << "\nVerify comparison operations\n";
m = pattern;
cout << " (m=pattern)>0\n";
assure( m > 0, "Unexpectedly, v is not greater than 0");
cout << " (m=pattern)>-pattern\n";
assert( m > -pattern );
cout << " (m=-pattern)<-pattern/2\n";
m -= 2*pattern;
assert( m < -pattern/2 );
cout << "\nAssign 2*pattern to m by repeating additions\n";
m = 0; m += pattern; m += pattern;
cout << "Assign 2*pattern to m1 by multiplying by two \n";
m1 = pattern; m1 *= 2;
verify_element_value(m1,2*pattern);
assert( m == m1 );
cout << "Multiply m1 by one half returning it to the 1*pattern\n";
m1 *= 1/2.;
verify_element_value(m1,pattern);
cout << "\nAssign -pattern to m and m1\n";
m.clear(); m -= pattern; m1 = -pattern;
verify_element_value(m,-pattern);
assert( m == m1 );
cout << "m = sqrt(sqr(m)); m1 = abs(m1); Now m and m1 have to be the same\n";
m.sqr();
verify_element_value(m,pattern*pattern);
m.sqrt();
verify_element_value(m,pattern);
m1.abs();
verify_element_value(m1,pattern);
assert( m == m1 );
cout << "\nCheck out to see that sin^2(x) + cos^2(x) = 1\n";
for(i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
m(i,j) = 4*PI/m.q_no_elems() * (i*m.q_ncols()+j);
m1 = m;
m.user_function(sin);
m1.user_function(cos);
m.sqr();
m1.sqr();
m += m1;
verify_element_value(m,1);
cout << "\nDone\n\n";
}
/*
*------------------------------------------------------------------------
* Test binary matrix element-by-element operations
*/
static void test_binary_ebe_op(const int rsize, const int csize)
{
const double pattern = 4.25;
register int i,j;
cout << "\n---> Test Binary Matrix element-by-element operations\n";
Matrix m(2,rsize+1,0,csize-1);
Matrix m1(m);
Matrix mp(m);
for(i=mp.q_row_lwb(); i<=mp.q_row_upb(); i++)
for(j=mp.q_col_lwb(); j<=mp.q_col_upb(); j++)
mp(i,j) = (i-m.q_nrows()/2.)*j*pattern;
cout << "\nVerify assignment of a matrix to the matrix\n";
m = pattern;
m1.clear();
m1 = m;
verify_element_value(m1,pattern);
assert( m1 == m );
cout << "\nAdding the matrix to itself, uniform pattern " << pattern << "\n";
m.clear(); m = pattern;
m1 = m; m1 += m1;
verify_element_value(m1,2*pattern);
cout << " subtracting two matrices ...\n";
m1 -= m;
verify_element_value(m1,pattern);
cout << " subtracting the matrix from itself\n";
m1 -= m1;
verify_element_value(m1,0);
cout << " adding two matrices together\n";
m1 += m;
verify_element_value(m1,pattern);
cout << "\nArithmetic operations on matrices with not the same elements\n";
cout << " adding mp to the zero matrix...\n";
m.clear(); m += mp;
verify_matrix_identity(m,mp);
m1 = m;
cout << " making m = 3*mp and m1 = 3*mp, via add() and succesive mult\n";
add(m,2,mp);
m1 += m1; m1 += mp;
verify_matrix_identity(m,m1);
cout << " clear both m and m1, by subtracting from itself and via add()\n";
m1 -= m1;
add(m,-3,mp);
verify_matrix_identity(m,m1);
cout << "\nTesting element-by-element multiplications and divisions\n";
cout << " squaring each element with sqr() and via multiplication\n";
m = mp; m1 = mp;
m.sqr();
element_mult(m1,m1);
verify_matrix_identity(m,m1);
cout << " compare (m = pattern^2)/pattern with pattern\n";
m = pattern; m1 = pattern;
m.sqr();
element_div(m,m1);
verify_element_value(m,pattern);
compare(m1,m,"Original vector and vector after squaring and dividing");
cout << "\nDone\n";
}
/*
*------------------------------------------------------------------------
* Verify the determinant evaluation
*/
static void test_determinant(const int msize)
{
cout << "\n---> Verify determinant evaluation\n"
"for a square matrix of size " << msize << "\n";
Matrix m(msize,msize);
cout << "\nCheck to see that the determinant of the unit matrix is one";
m.unit_matrix();
cout << "\n determinant is " << m.determinant();
assert( m.determinant() == 1 );
const double pattern = 2.5;
cout << "\nCheck the determinant for the matrix with " << pattern <<
"\n at the diagonal";
register int i,j;
for(i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
m(i,j) = ( i==j ? pattern : 0 );
cout << "\n determinant is " << m.determinant();
assert( m.determinant() == pow(pattern,m.q_nrows()) );
cout << "\nCheck the determinant for the matrix with " << pattern <<
"\n at the anti-diagonal";
for(i=m.q_row_lwb(); i<=m.q_row_upb(); i++)
for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
m(i,j) = ( i==(m.q_col_upb()+m.q_col_lwb()-j) ? pattern : 0 );
assert( m.determinant() == pow(pattern,m.q_nrows()) *
( m.q_nrows()*(m.q_nrows()-1)/2 & 1 ? -1 : 1 ) );
cout << "\nCheck the determinant for the singular matrix"
"\n defined as above with zero first row";
m.clear();
for(i=m.q_row_lwb()+1; i<=m.q_row_upb(); i++)
for(j=m.q_col_lwb(); j<=m.q_col_upb(); j++)
m(i,j) = ( i==(m.q_col_upb()+m.q_col_lwb()-j) ? pattern : 0 );
cout << "\n determinant is " << m.determinant();
assert( m.determinant() == 0 );
cout << "\nCheck out the determinant of the Hilbert matrix";
Matrix H(3,3);
H.hilbert_matrix();
cout << "\n 3x3 Hilbert matrix: exact determinant 1/2160 ";
cout << "\n computed 1/"<< 1/H.determinant();
H.resize_to(4,4);
H.hilbert_matrix();
cout << "\n 4x4 Hilbert matrix: exact determinant 1/6048000 ";
cout << "\n computed 1/"<< 1/H.determinant();
H.resize_to(5,5);
H.hilbert_matrix();
cout << "\n 5x5 Hilbert matrix: exact determinant 3.749295e-12";
cout << "\n computed "<< H.determinant();
H.resize_to(7,7);
H.hilbert_matrix();
cout << "\n 7x7 Hilbert matrix: exact determinant 4.8358e-25";
cout << "\n computed "<< H.determinant();
H.resize_to(9,9);
H.hilbert_matrix();
cout << "\n 9x9 Hilbert matrix: exact determinant 9.72023e-43";
cout << "\n computed "<< H.determinant();
H.resize_to(10,10);
H.hilbert_matrix();
cout << "\n 10x10 Hilbert matrix: exact determinant 2.16418e-53";
cout << "\n computed "<< H.determinant();
cout << "\nDone\n";
}
/*
*------------------------------------------------------------------------
* Verify the norm calculation
*/
static void test_norms(const int rsize, const int csize)
{
cout << "\n---> Verify norm calculations\n";
const double pattern = 10.25;
if( rsize % 2 == 1 || csize %2 == 1 )
_error("Sorry, size of the matrix to test must be even for this test\n");
Matrix m(rsize,csize);
Matrix m1(m);
cout << "\nAssign " << pattern << " to all the elements and check norms\n";
m = pattern;
cout << " 1. (col) norm should be pattern*nrows\n";
assert( m.norm_1() == pattern*m.q_nrows() );
assert( m.norm_1() == m.col_norm() );
cout << " Inf (row) norm should be pattern*ncols\n";
assert( m.norm_inf() == pattern*m.q_ncols() );
assert( m.norm_inf() == m.row_norm() );
cout << " Square of the Eucl norm has got to be pattern^2 * no_elems\n";
assert( m.e2_norm() == sqr(pattern)*m.q_no_elems() );
m1.clear();
assert( m.e2_norm() == e2_norm(m,m1) );
#if 0
double ap_step = 1;
double ap_a0 = -pattern;
int n = m.q_no_elems();
cout << "\nAssign the arithm progression with 1. term " << ap_a0 <<
"\nand the difference " << ap_step << "\n";
register int i;
for(i=v.q_lwb(); i<=v.q_upb(); i++)
v(i) = (i-v.q_lwb())*ap_step + ap_a0;
int l = min(max((int)ceil(-ap_a0/ap_step),0),n);
double norm = (2*ap_a0 + (l+n-1)*ap_step)/2*(n-l) +
(-2*ap_a0-(l-1)*ap_step)/2*l;
cout << " 1. norm should be " << norm << "\n";
assert( v.norm_1() == norm );
norm = n*( sqr(ap_a0) + ap_a0*ap_step*(n-1) + sqr(ap_step)*(n-1)*(2*n-1)/6);
cout << " Square of the 2. norm has got to be "
"n*[ a0^2 + a0*q*(n-1) + q^2/6*(n-1)*(2n-1) ], or " << norm << "\n";
assert( v.norm_2_sqr() == norm );
norm = max(abs(v(v.q_lwb())),abs(v(v.q_upb())));
cout << " Inf norm should be max(abs(a0),abs(a0+(n-1)*q)), ie " << norm <<
"\n";
assert( v.norm_inf() == norm );
m1.clear();
compare(m,m1,"Compare the matrix m with a zero matrix");
cout << "\nConstruct v1 to be orthogonal to v as v(n), -v(n-1), v(n-2)...\n";
for(i=0; i<v1.q_no_elems(); i++)
v1(i+v1.q_lwb()) = v(v.q_upb()-i) * ( i % 2 == 1 ? -1 : 1 );
cout << "||v1|| has got to be equal ||v|| regardless of the norm def\n";
assert( v1.norm_1() == v.norm_1() );
assert( v1.norm_2_sqr() == v.norm_2_sqr() );
assert( v1.norm_inf() == v.norm_inf() );
#endif
cout << "\nDone\n";
}
/*
*------------------------------------------------------------------------
* Root module
*/
main()
{
cout << "\n\n" << _Minuses <<
"\n\t\tVerify Operations on Matrices\n";
#if 1
test_allocation();
test_element_op(20,10);
test_binary_ebe_op(10,20);
test_norms(40,20);
#endif
test_determinant(20);
}